#include "custom_elements/stokes_3D.h"

namespace Kratos
{


void Stokes3D::ComputeGaussPointLHSContribution(BoundedMatrix<double,16,16>& lhs, const element_data<4,3>& data)
{
    const int nnodes = 4;
    const int dim = 3;

    const double rho = inner_prod(data.N, data.rho);
    const double& bdf0 = data.bdf0;

    //get constitutive matrix
    const Matrix& C = data.C;

    //get shape function values
    const BoundedMatrix<double,nnodes,dim>& DN = data.DN_DX;
    const array_1d<double,nnodes>& N = data.N;


    //compute an equivalent tau by Bitrans*c*Bi
    const double tau_denom =             (C(3,3) + C(4,4) + C(5,5))*(pow(DN(0,0), 2) + pow(DN(0,1), 2) + pow(DN(0,2), 2) + pow(DN(1,0), 2) + pow(DN(1,1), 2) + pow(DN(1,2), 2) + pow(DN(2,0), 2) + pow(DN(2,1), 2) + pow(DN(2,2), 2) + pow(DN(3,0), 2) + pow(DN(3,1), 2) + pow(DN(3,2), 2));

    const double tau1 = 1.0/(tau_denom);


    const double x0 =             bdf0*rho;
    const double x1 =             C(0,0)*DN(0,0) + C(0,3)*DN(0,1) + C(0,5)*DN(0,2);
    const double x2 =             C(0,3)*DN(0,0);
    const double x3 =             C(3,3)*DN(0,1) + C(3,5)*DN(0,2) + x2;
    const double x4 =             C(0,5)*DN(0,0);
    const double x5 =             C(3,5)*DN(0,1) + C(5,5)*DN(0,2) + x4;
    const double x6 =             C(0,1)*DN(0,1) + C(0,4)*DN(0,2) + x2;
    const double x7 =             C(1,3)*DN(0,1);
    const double x8 =             C(3,3)*DN(0,0) + C(3,4)*DN(0,2) + x7;
    const double x9 =             C(3,5)*DN(0,0);
    const double x10 =             C(4,5)*DN(0,2);
    const double x11 =             C(1,5)*DN(0,1) + x10 + x9;
    const double x12 =             C(0,2)*DN(0,2) + C(0,4)*DN(0,1) + x4;
    const double x13 =             C(3,4)*DN(0,1);
    const double x14 =             C(2,3)*DN(0,2) + x13 + x9;
    const double x15 =             C(2,5)*DN(0,2);
    const double x16 =             C(4,5)*DN(0,1) + C(5,5)*DN(0,0) + x15;
    const double x17 =             DN(0,0)*N[0];
    const double x18 =             C(0,0)*DN(1,0) + C(0,3)*DN(1,1) + C(0,5)*DN(1,2);
    const double x19 =             C(0,3)*DN(1,0);
    const double x20 =             C(3,3)*DN(1,1) + C(3,5)*DN(1,2) + x19;
    const double x21 =             C(0,5)*DN(1,0);
    const double x22 =             C(3,5)*DN(1,1) + C(5,5)*DN(1,2) + x21;
    const double x23 =             C(0,1)*DN(1,1) + C(0,4)*DN(1,2) + x19;
    const double x24 =             C(1,3)*DN(1,1);
    const double x25 =             C(3,3)*DN(1,0) + C(3,4)*DN(1,2) + x24;
    const double x26 =             C(3,5)*DN(1,0);
    const double x27 =             C(4,5)*DN(1,2);
    const double x28 =             C(1,5)*DN(1,1) + x26 + x27;
    const double x29 =             C(0,2)*DN(1,2) + C(0,4)*DN(1,1) + x21;
    const double x30 =             C(3,4)*DN(1,1);
    const double x31 =             C(2,3)*DN(1,2) + x26 + x30;
    const double x32 =             C(2,5)*DN(1,2);
    const double x33 =             C(4,5)*DN(1,1) + C(5,5)*DN(1,0) + x32;
    const double x34 =             DN(0,0)*N[1];
    const double x35 =             C(0,0)*DN(2,0) + C(0,3)*DN(2,1) + C(0,5)*DN(2,2);
    const double x36 =             C(0,3)*DN(2,0);
    const double x37 =             C(3,3)*DN(2,1) + C(3,5)*DN(2,2) + x36;
    const double x38 =             C(0,5)*DN(2,0);
    const double x39 =             C(3,5)*DN(2,1) + C(5,5)*DN(2,2) + x38;
    const double x40 =             C(0,1)*DN(2,1) + C(0,4)*DN(2,2) + x36;
    const double x41 =             C(1,3)*DN(2,1);
    const double x42 =             C(3,3)*DN(2,0) + C(3,4)*DN(2,2) + x41;
    const double x43 =             C(3,5)*DN(2,0);
    const double x44 =             C(4,5)*DN(2,2);
    const double x45 =             C(1,5)*DN(2,1) + x43 + x44;
    const double x46 =             C(0,2)*DN(2,2) + C(0,4)*DN(2,1) + x38;
    const double x47 =             C(3,4)*DN(2,1);
    const double x48 =             C(2,3)*DN(2,2) + x43 + x47;
    const double x49 =             C(2,5)*DN(2,2);
    const double x50 =             C(4,5)*DN(2,1) + C(5,5)*DN(2,0) + x49;
    const double x51 =             DN(0,0)*N[2];
    const double x52 =             C(0,0)*DN(3,0) + C(0,3)*DN(3,1) + C(0,5)*DN(3,2);
    const double x53 =             C(0,3)*DN(3,0);
    const double x54 =             C(3,3)*DN(3,1) + C(3,5)*DN(3,2) + x53;
    const double x55 =             C(0,5)*DN(3,0);
    const double x56 =             C(3,5)*DN(3,1) + C(5,5)*DN(3,2) + x55;
    const double x57 =             C(0,1)*DN(3,1) + C(0,4)*DN(3,2) + x53;
    const double x58 =             C(1,3)*DN(3,1);
    const double x59 =             C(3,3)*DN(3,0) + C(3,4)*DN(3,2) + x58;
    const double x60 =             C(3,5)*DN(3,0);
    const double x61 =             C(4,5)*DN(3,2);
    const double x62 =             C(1,5)*DN(3,1) + x60 + x61;
    const double x63 =             C(0,2)*DN(3,2) + C(0,4)*DN(3,1) + x55;
    const double x64 =             C(3,4)*DN(3,1);
    const double x65 =             C(2,3)*DN(3,2) + x60 + x64;
    const double x66 =             C(2,5)*DN(3,2);
    const double x67 =             C(4,5)*DN(3,1) + C(5,5)*DN(3,0) + x66;
    const double x68 =             DN(0,0)*N[3];
    const double x69 =             C(0,1)*DN(0,0) + C(1,5)*DN(0,2) + x7;
    const double x70 =             C(0,4)*DN(0,0) + x10 + x13;
    const double x71 =             C(1,1)*DN(0,1) + C(1,3)*DN(0,0) + C(1,4)*DN(0,2);
    const double x72 =             C(1,4)*DN(0,1);
    const double x73 =             C(3,4)*DN(0,0) + C(4,4)*DN(0,2) + x72;
    const double x74 =             C(1,2)*DN(0,2) + C(1,5)*DN(0,0) + x72;
    const double x75 =             C(2,4)*DN(0,2);
    const double x76 =             C(4,4)*DN(0,1) + C(4,5)*DN(0,0) + x75;
    const double x77 =             DN(0,1)*N[0];
    const double x78 =             C(0,1)*DN(1,0) + C(1,5)*DN(1,2) + x24;
    const double x79 =             C(0,4)*DN(1,0) + x27 + x30;
    const double x80 =             C(1,1)*DN(1,1) + C(1,3)*DN(1,0) + C(1,4)*DN(1,2);
    const double x81 =             C(1,4)*DN(1,1);
    const double x82 =             C(3,4)*DN(1,0) + C(4,4)*DN(1,2) + x81;
    const double x83 =             C(1,2)*DN(1,2) + C(1,5)*DN(1,0) + x81;
    const double x84 =             C(2,4)*DN(1,2);
    const double x85 =             C(4,4)*DN(1,1) + C(4,5)*DN(1,0) + x84;
    const double x86 =             DN(0,1)*N[1];
    const double x87 =             C(0,1)*DN(2,0) + C(1,5)*DN(2,2) + x41;
    const double x88 =             C(0,4)*DN(2,0) + x44 + x47;
    const double x89 =             C(1,1)*DN(2,1) + C(1,3)*DN(2,0) + C(1,4)*DN(2,2);
    const double x90 =             C(1,4)*DN(2,1);
    const double x91 =             C(3,4)*DN(2,0) + C(4,4)*DN(2,2) + x90;
    const double x92 =             C(1,2)*DN(2,2) + C(1,5)*DN(2,0) + x90;
    const double x93 =             C(2,4)*DN(2,2);
    const double x94 =             C(4,4)*DN(2,1) + C(4,5)*DN(2,0) + x93;
    const double x95 =             DN(0,1)*N[2];
    const double x96 =             C(0,1)*DN(3,0) + C(1,5)*DN(3,2) + x58;
    const double x97 =             C(0,4)*DN(3,0) + x61 + x64;
    const double x98 =             C(1,1)*DN(3,1) + C(1,3)*DN(3,0) + C(1,4)*DN(3,2);
    const double x99 =             C(1,4)*DN(3,1);
    const double x100 =             C(3,4)*DN(3,0) + C(4,4)*DN(3,2) + x99;
    const double x101 =             C(1,2)*DN(3,2) + C(1,5)*DN(3,0) + x99;
    const double x102 =             C(2,4)*DN(3,2);
    const double x103 =             C(4,4)*DN(3,1) + C(4,5)*DN(3,0) + x102;
    const double x104 =             DN(0,1)*N[3];
    const double x105 =             C(0,2)*DN(0,0) + C(2,3)*DN(0,1) + x15;
    const double x106 =             C(1,2)*DN(0,1) + C(2,3)*DN(0,0) + x75;
    const double x107 =             C(2,2)*DN(0,2) + C(2,4)*DN(0,1) + C(2,5)*DN(0,0);
    const double x108 =             DN(0,2)*N[0];
    const double x109 =             C(0,2)*DN(1,0) + C(2,3)*DN(1,1) + x32;
    const double x110 =             C(1,2)*DN(1,1) + C(2,3)*DN(1,0) + x84;
    const double x111 =             C(2,2)*DN(1,2) + C(2,4)*DN(1,1) + C(2,5)*DN(1,0);
    const double x112 =             DN(0,2)*N[1];
    const double x113 =             C(0,2)*DN(2,0) + C(2,3)*DN(2,1) + x49;
    const double x114 =             C(1,2)*DN(2,1) + C(2,3)*DN(2,0) + x93;
    const double x115 =             C(2,2)*DN(2,2) + C(2,4)*DN(2,1) + C(2,5)*DN(2,0);
    const double x116 =             DN(0,2)*N[2];
    const double x117 =             C(0,2)*DN(3,0) + C(2,3)*DN(3,1) + x66;
    const double x118 =             C(1,2)*DN(3,1) + C(2,3)*DN(3,0) + x102;
    const double x119 =             C(2,2)*DN(3,2) + C(2,4)*DN(3,1) + C(2,5)*DN(3,0);
    const double x120 =             DN(0,2)*N[3];
    const double x121 =             tau1*x0;
    const double x122 =             x121 + 1;
    const double x123 =             DN(1,0)*N[0];
    const double x124 =             DN(1,1)*N[0];
    const double x125 =             DN(1,2)*N[0];
    const double x126 =             tau1*(DN(0,0)*DN(1,0) + DN(0,1)*DN(1,1) + DN(0,2)*DN(1,2));
    const double x127 =             DN(2,0)*N[0];
    const double x128 =             DN(2,1)*N[0];
    const double x129 =             DN(2,2)*N[0];
    const double x130 =             tau1*(DN(0,0)*DN(2,0) + DN(0,1)*DN(2,1) + DN(0,2)*DN(2,2));
    const double x131 =             DN(3,0)*N[0];
    const double x132 =             DN(3,1)*N[0];
    const double x133 =             DN(3,2)*N[0];
    const double x134 =             tau1*(DN(0,0)*DN(3,0) + DN(0,1)*DN(3,1) + DN(0,2)*DN(3,2));
    const double x135 =             DN(1,0)*N[1];
    const double x136 =             DN(1,0)*N[2];
    const double x137 =             DN(1,0)*N[3];
    const double x138 =             DN(1,1)*N[1];
    const double x139 =             DN(1,1)*N[2];
    const double x140 =             DN(1,1)*N[3];
    const double x141 =             DN(1,2)*N[1];
    const double x142 =             DN(1,2)*N[2];
    const double x143 =             DN(1,2)*N[3];
    const double x144 =             DN(2,0)*N[1];
    const double x145 =             DN(2,1)*N[1];
    const double x146 =             DN(2,2)*N[1];
    const double x147 =             tau1*(DN(1,0)*DN(2,0) + DN(1,1)*DN(2,1) + DN(1,2)*DN(2,2));
    const double x148 =             DN(3,0)*N[1];
    const double x149 =             DN(3,1)*N[1];
    const double x150 =             DN(3,2)*N[1];
    const double x151 =             tau1*(DN(1,0)*DN(3,0) + DN(1,1)*DN(3,1) + DN(1,2)*DN(3,2));
    const double x152 =             DN(2,0)*N[2];
    const double x153 =             DN(2,0)*N[3];
    const double x154 =             DN(2,1)*N[2];
    const double x155 =             DN(2,1)*N[3];
    const double x156 =             DN(2,2)*N[2];
    const double x157 =             DN(2,2)*N[3];
    const double x158 =             DN(3,0)*N[2];
    const double x159 =             DN(3,1)*N[2];
    const double x160 =             DN(3,2)*N[2];
    const double x161 =             tau1*(DN(2,0)*DN(3,0) + DN(2,1)*DN(3,1) + DN(2,2)*DN(3,2));
    const double x162 =             DN(3,0)*N[3];
    const double x163 =             DN(3,1)*N[3];
    const double x164 =             DN(3,2)*N[3];
    lhs(0,0)=DN(0,0)*x1 + DN(0,1)*x3 + DN(0,2)*x5 + x0;
    lhs(0,1)=DN(0,0)*x6 + DN(0,1)*x8 + DN(0,2)*x11;
    lhs(0,2)=DN(0,0)*x12 + DN(0,1)*x14 + DN(0,2)*x16;
    lhs(0,3)=-x17;
    lhs(0,4)=DN(0,0)*x18 + DN(0,1)*x20 + DN(0,2)*x22;
    lhs(0,5)=DN(0,0)*x23 + DN(0,1)*x25 + DN(0,2)*x28;
    lhs(0,6)=DN(0,0)*x29 + DN(0,1)*x31 + DN(0,2)*x33;
    lhs(0,7)=-x34;
    lhs(0,8)=DN(0,0)*x35 + DN(0,1)*x37 + DN(0,2)*x39;
    lhs(0,9)=DN(0,0)*x40 + DN(0,1)*x42 + DN(0,2)*x45;
    lhs(0,10)=DN(0,0)*x46 + DN(0,1)*x48 + DN(0,2)*x50;
    lhs(0,11)=-x51;
    lhs(0,12)=DN(0,0)*x52 + DN(0,1)*x54 + DN(0,2)*x56;
    lhs(0,13)=DN(0,0)*x57 + DN(0,1)*x59 + DN(0,2)*x62;
    lhs(0,14)=DN(0,0)*x63 + DN(0,1)*x65 + DN(0,2)*x67;
    lhs(0,15)=-x68;
    lhs(1,0)=DN(0,0)*x3 + DN(0,1)*x69 + DN(0,2)*x70;
    lhs(1,1)=DN(0,0)*x8 + DN(0,1)*x71 + DN(0,2)*x73 + x0;
    lhs(1,2)=DN(0,0)*x14 + DN(0,1)*x74 + DN(0,2)*x76;
    lhs(1,3)=-x77;
    lhs(1,4)=DN(0,0)*x20 + DN(0,1)*x78 + DN(0,2)*x79;
    lhs(1,5)=DN(0,0)*x25 + DN(0,1)*x80 + DN(0,2)*x82;
    lhs(1,6)=DN(0,0)*x31 + DN(0,1)*x83 + DN(0,2)*x85;
    lhs(1,7)=-x86;
    lhs(1,8)=DN(0,0)*x37 + DN(0,1)*x87 + DN(0,2)*x88;
    lhs(1,9)=DN(0,0)*x42 + DN(0,1)*x89 + DN(0,2)*x91;
    lhs(1,10)=DN(0,0)*x48 + DN(0,1)*x92 + DN(0,2)*x94;
    lhs(1,11)=-x95;
    lhs(1,12)=DN(0,0)*x54 + DN(0,1)*x96 + DN(0,2)*x97;
    lhs(1,13)=DN(0,0)*x59 + DN(0,1)*x98 + DN(0,2)*x100;
    lhs(1,14)=DN(0,0)*x65 + DN(0,1)*x101 + DN(0,2)*x103;
    lhs(1,15)=-x104;
    lhs(2,0)=DN(0,0)*x5 + DN(0,1)*x70 + DN(0,2)*x105;
    lhs(2,1)=DN(0,0)*x11 + DN(0,1)*x73 + DN(0,2)*x106;
    lhs(2,2)=DN(0,0)*x16 + DN(0,1)*x76 + DN(0,2)*x107 + x0;
    lhs(2,3)=-x108;
    lhs(2,4)=DN(0,0)*x22 + DN(0,1)*x79 + DN(0,2)*x109;
    lhs(2,5)=DN(0,0)*x28 + DN(0,1)*x82 + DN(0,2)*x110;
    lhs(2,6)=DN(0,0)*x33 + DN(0,1)*x85 + DN(0,2)*x111;
    lhs(2,7)=-x112;
    lhs(2,8)=DN(0,0)*x39 + DN(0,1)*x88 + DN(0,2)*x113;
    lhs(2,9)=DN(0,0)*x45 + DN(0,1)*x91 + DN(0,2)*x114;
    lhs(2,10)=DN(0,0)*x50 + DN(0,1)*x94 + DN(0,2)*x115;
    lhs(2,11)=-x116;
    lhs(2,12)=DN(0,0)*x56 + DN(0,1)*x97 + DN(0,2)*x117;
    lhs(2,13)=DN(0,0)*x62 + DN(0,1)*x100 + DN(0,2)*x118;
    lhs(2,14)=DN(0,0)*x67 + DN(0,1)*x103 + DN(0,2)*x119;
    lhs(2,15)=-x120;
    lhs(3,0)=x122*x17;
    lhs(3,1)=x122*x77;
    lhs(3,2)=x108*x122;
    lhs(3,3)=tau1*(pow(DN(0,0), 2) + pow(DN(0,1), 2) + pow(DN(0,2), 2));
    lhs(3,4)=x121*x34 + x123;
    lhs(3,5)=x121*x86 + x124;
    lhs(3,6)=x112*x121 + x125;
    lhs(3,7)=x126;
    lhs(3,8)=x121*x51 + x127;
    lhs(3,9)=x121*x95 + x128;
    lhs(3,10)=x116*x121 + x129;
    lhs(3,11)=x130;
    lhs(3,12)=x121*x68 + x131;
    lhs(3,13)=x104*x121 + x132;
    lhs(3,14)=x120*x121 + x133;
    lhs(3,15)=x134;
    lhs(4,0)=DN(1,0)*x1 + DN(1,1)*x3 + DN(1,2)*x5;
    lhs(4,1)=DN(1,0)*x6 + DN(1,1)*x8 + DN(1,2)*x11;
    lhs(4,2)=DN(1,0)*x12 + DN(1,1)*x14 + DN(1,2)*x16;
    lhs(4,3)=-x123;
    lhs(4,4)=DN(1,0)*x18 + DN(1,1)*x20 + DN(1,2)*x22 + x0;
    lhs(4,5)=DN(1,0)*x23 + DN(1,1)*x25 + DN(1,2)*x28;
    lhs(4,6)=DN(1,0)*x29 + DN(1,1)*x31 + DN(1,2)*x33;
    lhs(4,7)=-x135;
    lhs(4,8)=DN(1,0)*x35 + DN(1,1)*x37 + DN(1,2)*x39;
    lhs(4,9)=DN(1,0)*x40 + DN(1,1)*x42 + DN(1,2)*x45;
    lhs(4,10)=DN(1,0)*x46 + DN(1,1)*x48 + DN(1,2)*x50;
    lhs(4,11)=-x136;
    lhs(4,12)=DN(1,0)*x52 + DN(1,1)*x54 + DN(1,2)*x56;
    lhs(4,13)=DN(1,0)*x57 + DN(1,1)*x59 + DN(1,2)*x62;
    lhs(4,14)=DN(1,0)*x63 + DN(1,1)*x65 + DN(1,2)*x67;
    lhs(4,15)=-x137;
    lhs(5,0)=DN(1,0)*x3 + DN(1,1)*x69 + DN(1,2)*x70;
    lhs(5,1)=DN(1,0)*x8 + DN(1,1)*x71 + DN(1,2)*x73;
    lhs(5,2)=DN(1,0)*x14 + DN(1,1)*x74 + DN(1,2)*x76;
    lhs(5,3)=-x124;
    lhs(5,4)=DN(1,0)*x20 + DN(1,1)*x78 + DN(1,2)*x79;
    lhs(5,5)=DN(1,0)*x25 + DN(1,1)*x80 + DN(1,2)*x82 + x0;
    lhs(5,6)=DN(1,0)*x31 + DN(1,1)*x83 + DN(1,2)*x85;
    lhs(5,7)=-x138;
    lhs(5,8)=DN(1,0)*x37 + DN(1,1)*x87 + DN(1,2)*x88;
    lhs(5,9)=DN(1,0)*x42 + DN(1,1)*x89 + DN(1,2)*x91;
    lhs(5,10)=DN(1,0)*x48 + DN(1,1)*x92 + DN(1,2)*x94;
    lhs(5,11)=-x139;
    lhs(5,12)=DN(1,0)*x54 + DN(1,1)*x96 + DN(1,2)*x97;
    lhs(5,13)=DN(1,0)*x59 + DN(1,1)*x98 + DN(1,2)*x100;
    lhs(5,14)=DN(1,0)*x65 + DN(1,1)*x101 + DN(1,2)*x103;
    lhs(5,15)=-x140;
    lhs(6,0)=DN(1,0)*x5 + DN(1,1)*x70 + DN(1,2)*x105;
    lhs(6,1)=DN(1,0)*x11 + DN(1,1)*x73 + DN(1,2)*x106;
    lhs(6,2)=DN(1,0)*x16 + DN(1,1)*x76 + DN(1,2)*x107;
    lhs(6,3)=-x125;
    lhs(6,4)=DN(1,0)*x22 + DN(1,1)*x79 + DN(1,2)*x109;
    lhs(6,5)=DN(1,0)*x28 + DN(1,1)*x82 + DN(1,2)*x110;
    lhs(6,6)=DN(1,0)*x33 + DN(1,1)*x85 + DN(1,2)*x111 + x0;
    lhs(6,7)=-x141;
    lhs(6,8)=DN(1,0)*x39 + DN(1,1)*x88 + DN(1,2)*x113;
    lhs(6,9)=DN(1,0)*x45 + DN(1,1)*x91 + DN(1,2)*x114;
    lhs(6,10)=DN(1,0)*x50 + DN(1,1)*x94 + DN(1,2)*x115;
    lhs(6,11)=-x142;
    lhs(6,12)=DN(1,0)*x56 + DN(1,1)*x97 + DN(1,2)*x117;
    lhs(6,13)=DN(1,0)*x62 + DN(1,1)*x100 + DN(1,2)*x118;
    lhs(6,14)=DN(1,0)*x67 + DN(1,1)*x103 + DN(1,2)*x119;
    lhs(6,15)=-x143;
    lhs(7,0)=x121*x123 + x34;
    lhs(7,1)=x121*x124 + x86;
    lhs(7,2)=x112 + x121*x125;
    lhs(7,3)=x126;
    lhs(7,4)=x122*x135;
    lhs(7,5)=x122*x138;
    lhs(7,6)=x122*x141;
    lhs(7,7)=tau1*(pow(DN(1,0), 2) + pow(DN(1,1), 2) + pow(DN(1,2), 2));
    lhs(7,8)=x121*x136 + x144;
    lhs(7,9)=x121*x139 + x145;
    lhs(7,10)=x121*x142 + x146;
    lhs(7,11)=x147;
    lhs(7,12)=x121*x137 + x148;
    lhs(7,13)=x121*x140 + x149;
    lhs(7,14)=x121*x143 + x150;
    lhs(7,15)=x151;
    lhs(8,0)=DN(2,0)*x1 + DN(2,1)*x3 + DN(2,2)*x5;
    lhs(8,1)=DN(2,0)*x6 + DN(2,1)*x8 + DN(2,2)*x11;
    lhs(8,2)=DN(2,0)*x12 + DN(2,1)*x14 + DN(2,2)*x16;
    lhs(8,3)=-x127;
    lhs(8,4)=DN(2,0)*x18 + DN(2,1)*x20 + DN(2,2)*x22;
    lhs(8,5)=DN(2,0)*x23 + DN(2,1)*x25 + DN(2,2)*x28;
    lhs(8,6)=DN(2,0)*x29 + DN(2,1)*x31 + DN(2,2)*x33;
    lhs(8,7)=-x144;
    lhs(8,8)=DN(2,0)*x35 + DN(2,1)*x37 + DN(2,2)*x39 + x0;
    lhs(8,9)=DN(2,0)*x40 + DN(2,1)*x42 + DN(2,2)*x45;
    lhs(8,10)=DN(2,0)*x46 + DN(2,1)*x48 + DN(2,2)*x50;
    lhs(8,11)=-x152;
    lhs(8,12)=DN(2,0)*x52 + DN(2,1)*x54 + DN(2,2)*x56;
    lhs(8,13)=DN(2,0)*x57 + DN(2,1)*x59 + DN(2,2)*x62;
    lhs(8,14)=DN(2,0)*x63 + DN(2,1)*x65 + DN(2,2)*x67;
    lhs(8,15)=-x153;
    lhs(9,0)=DN(2,0)*x3 + DN(2,1)*x69 + DN(2,2)*x70;
    lhs(9,1)=DN(2,0)*x8 + DN(2,1)*x71 + DN(2,2)*x73;
    lhs(9,2)=DN(2,0)*x14 + DN(2,1)*x74 + DN(2,2)*x76;
    lhs(9,3)=-x128;
    lhs(9,4)=DN(2,0)*x20 + DN(2,1)*x78 + DN(2,2)*x79;
    lhs(9,5)=DN(2,0)*x25 + DN(2,1)*x80 + DN(2,2)*x82;
    lhs(9,6)=DN(2,0)*x31 + DN(2,1)*x83 + DN(2,2)*x85;
    lhs(9,7)=-x145;
    lhs(9,8)=DN(2,0)*x37 + DN(2,1)*x87 + DN(2,2)*x88;
    lhs(9,9)=DN(2,0)*x42 + DN(2,1)*x89 + DN(2,2)*x91 + x0;
    lhs(9,10)=DN(2,0)*x48 + DN(2,1)*x92 + DN(2,2)*x94;
    lhs(9,11)=-x154;
    lhs(9,12)=DN(2,0)*x54 + DN(2,1)*x96 + DN(2,2)*x97;
    lhs(9,13)=DN(2,0)*x59 + DN(2,1)*x98 + DN(2,2)*x100;
    lhs(9,14)=DN(2,0)*x65 + DN(2,1)*x101 + DN(2,2)*x103;
    lhs(9,15)=-x155;
    lhs(10,0)=DN(2,0)*x5 + DN(2,1)*x70 + DN(2,2)*x105;
    lhs(10,1)=DN(2,0)*x11 + DN(2,1)*x73 + DN(2,2)*x106;
    lhs(10,2)=DN(2,0)*x16 + DN(2,1)*x76 + DN(2,2)*x107;
    lhs(10,3)=-x129;
    lhs(10,4)=DN(2,0)*x22 + DN(2,1)*x79 + DN(2,2)*x109;
    lhs(10,5)=DN(2,0)*x28 + DN(2,1)*x82 + DN(2,2)*x110;
    lhs(10,6)=DN(2,0)*x33 + DN(2,1)*x85 + DN(2,2)*x111;
    lhs(10,7)=-x146;
    lhs(10,8)=DN(2,0)*x39 + DN(2,1)*x88 + DN(2,2)*x113;
    lhs(10,9)=DN(2,0)*x45 + DN(2,1)*x91 + DN(2,2)*x114;
    lhs(10,10)=DN(2,0)*x50 + DN(2,1)*x94 + DN(2,2)*x115 + x0;
    lhs(10,11)=-x156;
    lhs(10,12)=DN(2,0)*x56 + DN(2,1)*x97 + DN(2,2)*x117;
    lhs(10,13)=DN(2,0)*x62 + DN(2,1)*x100 + DN(2,2)*x118;
    lhs(10,14)=DN(2,0)*x67 + DN(2,1)*x103 + DN(2,2)*x119;
    lhs(10,15)=-x157;
    lhs(11,0)=x121*x127 + x51;
    lhs(11,1)=x121*x128 + x95;
    lhs(11,2)=x116 + x121*x129;
    lhs(11,3)=x130;
    lhs(11,4)=x121*x144 + x136;
    lhs(11,5)=x121*x145 + x139;
    lhs(11,6)=x121*x146 + x142;
    lhs(11,7)=x147;
    lhs(11,8)=x122*x152;
    lhs(11,9)=x122*x154;
    lhs(11,10)=x122*x156;
    lhs(11,11)=tau1*(pow(DN(2,0), 2) + pow(DN(2,1), 2) + pow(DN(2,2), 2));
    lhs(11,12)=x121*x153 + x158;
    lhs(11,13)=x121*x155 + x159;
    lhs(11,14)=x121*x157 + x160;
    lhs(11,15)=x161;
    lhs(12,0)=DN(3,0)*x1 + DN(3,1)*x3 + DN(3,2)*x5;
    lhs(12,1)=DN(3,0)*x6 + DN(3,1)*x8 + DN(3,2)*x11;
    lhs(12,2)=DN(3,0)*x12 + DN(3,1)*x14 + DN(3,2)*x16;
    lhs(12,3)=-x131;
    lhs(12,4)=DN(3,0)*x18 + DN(3,1)*x20 + DN(3,2)*x22;
    lhs(12,5)=DN(3,0)*x23 + DN(3,1)*x25 + DN(3,2)*x28;
    lhs(12,6)=DN(3,0)*x29 + DN(3,1)*x31 + DN(3,2)*x33;
    lhs(12,7)=-x148;
    lhs(12,8)=DN(3,0)*x35 + DN(3,1)*x37 + DN(3,2)*x39;
    lhs(12,9)=DN(3,0)*x40 + DN(3,1)*x42 + DN(3,2)*x45;
    lhs(12,10)=DN(3,0)*x46 + DN(3,1)*x48 + DN(3,2)*x50;
    lhs(12,11)=-x158;
    lhs(12,12)=DN(3,0)*x52 + DN(3,1)*x54 + DN(3,2)*x56 + x0;
    lhs(12,13)=DN(3,0)*x57 + DN(3,1)*x59 + DN(3,2)*x62;
    lhs(12,14)=DN(3,0)*x63 + DN(3,1)*x65 + DN(3,2)*x67;
    lhs(12,15)=-x162;
    lhs(13,0)=DN(3,0)*x3 + DN(3,1)*x69 + DN(3,2)*x70;
    lhs(13,1)=DN(3,0)*x8 + DN(3,1)*x71 + DN(3,2)*x73;
    lhs(13,2)=DN(3,0)*x14 + DN(3,1)*x74 + DN(3,2)*x76;
    lhs(13,3)=-x132;
    lhs(13,4)=DN(3,0)*x20 + DN(3,1)*x78 + DN(3,2)*x79;
    lhs(13,5)=DN(3,0)*x25 + DN(3,1)*x80 + DN(3,2)*x82;
    lhs(13,6)=DN(3,0)*x31 + DN(3,1)*x83 + DN(3,2)*x85;
    lhs(13,7)=-x149;
    lhs(13,8)=DN(3,0)*x37 + DN(3,1)*x87 + DN(3,2)*x88;
    lhs(13,9)=DN(3,0)*x42 + DN(3,1)*x89 + DN(3,2)*x91;
    lhs(13,10)=DN(3,0)*x48 + DN(3,1)*x92 + DN(3,2)*x94;
    lhs(13,11)=-x159;
    lhs(13,12)=DN(3,0)*x54 + DN(3,1)*x96 + DN(3,2)*x97;
    lhs(13,13)=DN(3,0)*x59 + DN(3,1)*x98 + DN(3,2)*x100 + x0;
    lhs(13,14)=DN(3,0)*x65 + DN(3,1)*x101 + DN(3,2)*x103;
    lhs(13,15)=-x163;
    lhs(14,0)=DN(3,0)*x5 + DN(3,1)*x70 + DN(3,2)*x105;
    lhs(14,1)=DN(3,0)*x11 + DN(3,1)*x73 + DN(3,2)*x106;
    lhs(14,2)=DN(3,0)*x16 + DN(3,1)*x76 + DN(3,2)*x107;
    lhs(14,3)=-x133;
    lhs(14,4)=DN(3,0)*x22 + DN(3,1)*x79 + DN(3,2)*x109;
    lhs(14,5)=DN(3,0)*x28 + DN(3,1)*x82 + DN(3,2)*x110;
    lhs(14,6)=DN(3,0)*x33 + DN(3,1)*x85 + DN(3,2)*x111;
    lhs(14,7)=-x150;
    lhs(14,8)=DN(3,0)*x39 + DN(3,1)*x88 + DN(3,2)*x113;
    lhs(14,9)=DN(3,0)*x45 + DN(3,1)*x91 + DN(3,2)*x114;
    lhs(14,10)=DN(3,0)*x50 + DN(3,1)*x94 + DN(3,2)*x115;
    lhs(14,11)=-x160;
    lhs(14,12)=DN(3,0)*x56 + DN(3,1)*x97 + DN(3,2)*x117;
    lhs(14,13)=DN(3,0)*x62 + DN(3,1)*x100 + DN(3,2)*x118;
    lhs(14,14)=DN(3,0)*x67 + DN(3,1)*x103 + DN(3,2)*x119 + x0;
    lhs(14,15)=-x164;
    lhs(15,0)=x121*x131 + x68;
    lhs(15,1)=x104 + x121*x132;
    lhs(15,2)=x120 + x121*x133;
    lhs(15,3)=x134;
    lhs(15,4)=x121*x148 + x137;
    lhs(15,5)=x121*x149 + x140;
    lhs(15,6)=x121*x150 + x143;
    lhs(15,7)=x151;
    lhs(15,8)=x121*x158 + x153;
    lhs(15,9)=x121*x159 + x155;
    lhs(15,10)=x121*x160 + x157;
    lhs(15,11)=x161;
    lhs(15,12)=x122*x162;
    lhs(15,13)=x122*x163;
    lhs(15,14)=x122*x164;
    lhs(15,15)=tau1*(pow(DN(3,0), 2) + pow(DN(3,1), 2) + pow(DN(3,2), 2));


}

void Stokes3D::ComputeGaussPointRHSContribution(array_1d<double,16>& rhs, const element_data<4,3>& data)
{
    const int nnodes = 4;
    const int dim = 3;

    const double rho = inner_prod(data.N, data.rho);
    const double& bdf0 = data.bdf0;
    const double& bdf1 = data.bdf1;
    const double& bdf2 = data.bdf2;

    const BoundedMatrix<double,nnodes,dim>& v = data.v;
    const BoundedMatrix<double,nnodes,dim>& vn = data.vn;
    const BoundedMatrix<double,nnodes,dim>& vnn = data.vnn;
    const BoundedMatrix<double,nnodes,dim>& f = data.f;
    const array_1d<double,nnodes>& p = data.p;

    //get constitutive matrix
    const Matrix& C = data.C;
    const Vector& stress = data.stress;

    //get shape function values
    const BoundedMatrix<double,nnodes,dim>& DN = data.DN_DX;
    const array_1d<double,nnodes>& N = data.N;

    //compute an equivalent tau by Bitrans*c*Bi
    const double tau_denom =             (C(3,3) + C(4,4) + C(5,5))*(pow(DN(0,0), 2) + pow(DN(0,1), 2) + pow(DN(0,2), 2) + pow(DN(1,0), 2) + pow(DN(1,1), 2) + pow(DN(1,2), 2) + pow(DN(2,0), 2) + pow(DN(2,1), 2) + pow(DN(2,2), 2) + pow(DN(3,0), 2) + pow(DN(3,1), 2) + pow(DN(3,2), 2));

    const double tau1 = 1.0/(tau_denom);

    //auxiliary variables used in the calculation of the RHS
    const array_1d<double,dim> fgauss = prod(trans(f), N);
    const array_1d<double,dim> vgauss = prod(trans(v), N);
    const array_1d<double,dim> grad_p = prod(trans(DN), p);
    const double pgauss = inner_prod(N,p);

    array_1d<double,dim> acch = bdf0*vgauss;
    noalias(acch) += bdf1*prod(trans(vn), N);
    noalias(acch) += bdf2*prod(trans(vnn), N);

    rhs[0]=DN(0,0)*pgauss - DN(0,0)*stress[0] - DN(0,1)*stress[3] - DN(0,2)*stress[5] + N[0]*fgauss[0] - rho*(bdf0*v(0,0) + bdf1*vn(0,0) + bdf2*vnn(0,0));
    rhs[1]=-DN(0,0)*stress[3] + DN(0,1)*pgauss - DN(0,1)*stress[1] - DN(0,2)*stress[4] + N[0]*fgauss[1] - rho*(bdf0*v(0,1) + bdf1*vn(0,1) + bdf2*vnn(0,1));
    rhs[2]=-DN(0,0)*stress[5] - DN(0,1)*stress[4] + DN(0,2)*pgauss - DN(0,2)*stress[2] + N[0]*fgauss[2] - rho*(bdf0*v(0,2) + bdf1*vn(0,2) + bdf2*vnn(0,2));
    rhs[3]=-DN(0,0)*tau1*(acch[0]*rho - fgauss[0] + grad_p[0]) - DN(0,1)*tau1*(acch[1]*rho - fgauss[1] + grad_p[1]) - DN(0,2)*tau1*(acch[2]*rho - fgauss[2] + grad_p[2]) - N[0]*(DN(0,0)*v(0,0) + DN(0,1)*v(0,1) + DN(0,2)*v(0,2) + DN(1,0)*v(1,0) + DN(1,1)*v(1,1) + DN(1,2)*v(1,2) + DN(2,0)*v(2,0) + DN(2,1)*v(2,1) + DN(2,2)*v(2,2) + DN(3,0)*v(3,0) + DN(3,1)*v(3,1) + DN(3,2)*v(3,2));
    rhs[4]=DN(1,0)*pgauss - DN(1,0)*stress[0] - DN(1,1)*stress[3] - DN(1,2)*stress[5] + N[1]*fgauss[0] - rho*(bdf0*v(1,0) + bdf1*vn(1,0) + bdf2*vnn(1,0));
    rhs[5]=-DN(1,0)*stress[3] + DN(1,1)*pgauss - DN(1,1)*stress[1] - DN(1,2)*stress[4] + N[1]*fgauss[1] - rho*(bdf0*v(1,1) + bdf1*vn(1,1) + bdf2*vnn(1,1));
    rhs[6]=-DN(1,0)*stress[5] - DN(1,1)*stress[4] + DN(1,2)*pgauss - DN(1,2)*stress[2] + N[1]*fgauss[2] - rho*(bdf0*v(1,2) + bdf1*vn(1,2) + bdf2*vnn(1,2));
    rhs[7]=-DN(1,0)*tau1*(acch[0]*rho - fgauss[0] + grad_p[0]) - DN(1,1)*tau1*(acch[1]*rho - fgauss[1] + grad_p[1]) - DN(1,2)*tau1*(acch[2]*rho - fgauss[2] + grad_p[2]) - N[1]*(DN(0,0)*v(0,0) + DN(0,1)*v(0,1) + DN(0,2)*v(0,2) + DN(1,0)*v(1,0) + DN(1,1)*v(1,1) + DN(1,2)*v(1,2) + DN(2,0)*v(2,0) + DN(2,1)*v(2,1) + DN(2,2)*v(2,2) + DN(3,0)*v(3,0) + DN(3,1)*v(3,1) + DN(3,2)*v(3,2));
    rhs[8]=DN(2,0)*pgauss - DN(2,0)*stress[0] - DN(2,1)*stress[3] - DN(2,2)*stress[5] + N[2]*fgauss[0] - rho*(bdf0*v(2,0) + bdf1*vn(2,0) + bdf2*vnn(2,0));
    rhs[9]=-DN(2,0)*stress[3] + DN(2,1)*pgauss - DN(2,1)*stress[1] - DN(2,2)*stress[4] + N[2]*fgauss[1] - rho*(bdf0*v(2,1) + bdf1*vn(2,1) + bdf2*vnn(2,1));
    rhs[10]=-DN(2,0)*stress[5] - DN(2,1)*stress[4] + DN(2,2)*pgauss - DN(2,2)*stress[2] + N[2]*fgauss[2] - rho*(bdf0*v(2,2) + bdf1*vn(2,2) + bdf2*vnn(2,2));
    rhs[11]=-DN(2,0)*tau1*(acch[0]*rho - fgauss[0] + grad_p[0]) - DN(2,1)*tau1*(acch[1]*rho - fgauss[1] + grad_p[1]) - DN(2,2)*tau1*(acch[2]*rho - fgauss[2] + grad_p[2]) - N[2]*(DN(0,0)*v(0,0) + DN(0,1)*v(0,1) + DN(0,2)*v(0,2) + DN(1,0)*v(1,0) + DN(1,1)*v(1,1) + DN(1,2)*v(1,2) + DN(2,0)*v(2,0) + DN(2,1)*v(2,1) + DN(2,2)*v(2,2) + DN(3,0)*v(3,0) + DN(3,1)*v(3,1) + DN(3,2)*v(3,2));
    rhs[12]=DN(3,0)*pgauss - DN(3,0)*stress[0] - DN(3,1)*stress[3] - DN(3,2)*stress[5] + N[3]*fgauss[0] - rho*(bdf0*v(3,0) + bdf1*vn(3,0) + bdf2*vnn(3,0));
    rhs[13]=-DN(3,0)*stress[3] + DN(3,1)*pgauss - DN(3,1)*stress[1] - DN(3,2)*stress[4] + N[3]*fgauss[1] - rho*(bdf0*v(3,1) + bdf1*vn(3,1) + bdf2*vnn(3,1));
    rhs[14]=-DN(3,0)*stress[5] - DN(3,1)*stress[4] + DN(3,2)*pgauss - DN(3,2)*stress[2] + N[3]*fgauss[2] - rho*(bdf0*v(3,2) + bdf1*vn(3,2) + bdf2*vnn(3,2));
    rhs[15]=-DN(3,0)*tau1*(acch[0]*rho - fgauss[0] + grad_p[0]) - DN(3,1)*tau1*(acch[1]*rho - fgauss[1] + grad_p[1]) - DN(3,2)*tau1*(acch[2]*rho - fgauss[2] + grad_p[2]) - N[3]*(DN(0,0)*v(0,0) + DN(0,1)*v(0,1) + DN(0,2)*v(0,2) + DN(1,0)*v(1,0) + DN(1,1)*v(1,1) + DN(1,2)*v(1,2) + DN(2,0)*v(2,0) + DN(2,1)*v(2,1) + DN(2,2)*v(2,2) + DN(3,0)*v(3,0) + DN(3,1)*v(3,1) + DN(3,2)*v(3,2));

}

}

